######################### Load Libraries #######################################

library(foreign)

library(Zelig)

library(mvtnorm)

############################ Functions #########################################

# function that displays dimensions of a matrix

dims <- function (x) 
{
    if (is.vector(x)) 
        return(length(x))
    else return(dim(x))
}

############################ Data Setup ########################################

# set working directory
setwd("~/_data/")

# load LPR data set that includes variables from other data sets
load("lprdata_plus.RData")

################################### MS Table 2 #################################

model1.lm <- lm(p_new ~ lpr + lpr2 + cgdp_all + import_all + lnpop + growth_all + dxr_pwt + usainf_consumer + dec1980 + dec1990 + dec2000, data = lprprices[lprprices$irreg==0,])

model2.lm <- lm(p_new ~ vmarginall + cgdp_all + import_all + lnpop + growth_all + dxr_pwt + usainf_consumer + dec1980 + dec1990 + dec2000, data = lprprices[lprprices$irreg==0,])

model3.lm <- lm(p_new ~ seatsharegap + cgdp_all + import_all + lnpop + growth_all + dxr_pwt + usainf_consumer + dec1980 + dec1990 + dec2000, data = lprprices[lprprices$irreg==0,])

################################## Simulations #################################

set.seed(1234)

beta <- model1.lm$coefficients

VCV <- sandwich(model1.lm)

S <- rmvnorm(10000, beta, VCV)

x <- matrix(NA,161,dims(model.matrix(model1.lm))[2])

x.means <- apply(model.matrix(model1.lm)[,c(1,4:9)], 2, mean)

for (i in 1:161) {

	x[i,c(1,4:9)] <- x.means
	
}

x[,2] <- seq(0,0.8,0.005)

x[,3] <- x[,2]^2

x[,10] <- rep(0,161)

x[,11] <- rep(0,161)

x[,12] <- rep(1,161)


ev <- matrix(NA,10000,161)

for (i in 1:161) { 

	ev[,i] <- S%*%x[i,]
	
}

################################## MS Figure 2 #################################

# set working directory for graphs
setwd("~/_manuscript_tabfig/")

# save graph to pdf file
pdf("lpr_effect_on_prices.pdf")

par(fig=c(0,1,0.4,1), new=TRUE)
plot(1, type="n", xlim=c(0,.8), ylim=c(80,100), 
	xlab="Loss Probability", ylab="Estimated Values for Real Prices")
lines(seq(0,0.8,0.005), apply(ev,2, mean), col="blue")
lines(seq(0,0.8,0.005), apply(ev,2, quantile, prob=0.0825), col="blue", lty=2)
lines(seq(0,0.8,0.005), apply(ev,2, quantile, prob=0.9175), col="blue", lty=2)

par(fig=c(0,1,0,0.5), new=TRUE)
plot(density(model.matrix(model1.lm)[,2]),main="",xlab="Empirical Support for Loss Probability", col="red")

dev.off()